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We propose a method for reconstruction of the density ma- 
trix from measurable time-dependent (probability) distribu- 
tions of physical quantities. The applicability of the method 
based on least-squares inversion is - compared with other 
methods - very universal. It can be used to reconstruct 
quantum states of various systems, such as harmonic and 
and anharmonic oscillators including molecular vibrations in 
vibronic transitions and damped motion. It also enables 
one to take into account various specific features of experi- 
ments, such as limited sets of data and data smearing owing 
to limited resolution. To illustrate the method, we consider 
a Morse oscillator and give a comparison with other state- 
reconstruction methods suggested recently. 
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I. INTRODUCTION 

During the last years the problem of quantum-state 
measurement has been of increasing interest. Advances 
in the experimental techniques, which have allowed the 
physicists to measure not only single observables of a 
quantum-mechanical system but - for certain systems 
- the quantum state as a whole, have stimulated both 
experimental and theoretical research. A number of pro- 
posals for measuring quantum states have been made and 
various quantum-mechanical systems have been consid- 
ered. 

In quantum optics balanced homodyning has been a 
very fruitful method for tomographical reconstruction 
|jj of the quantum state of optical fields. Measuring 
the quadrature-component statistics and using inverse 
Radon transform techniques, the method was success- 
fully applied to the determination of the Wigner func- 
tion of single-mode signal fields The method (also 
called optical homodyne tomography) has also been ex- 
tended to direct sampling of the density matrix in a 
quadrature-component basis Q and in the Fock basis 
UK- Tomographical methods can also be used for re- 
constructing the quantum state of matter systems, such 
as molecular vibrations || or the transverse motion of 
an atom beam [Q. In the case of molecular vibrations 
the time- and frequency-resolved fluorescence spectrum 
of the molecule plays the role of the quadrature compo- 



nents, provided that the vibrational frequencies in the 
two electronic states are almost equal to each other and 
the vibrational motion is approximately harmonic. Fur- 
ther, tomographical methods have been considered in the 
reconstruction of the (harmonic) center-of-mass motion 
of trapped atoms Q . 

Including in the balanced detection scheme additional 
vacuum inputs, the quantum state of optical fields can be 
directly measured in terms of phase-space functions 
Using unbalanced homodyning, the displaced photon- 
number statistics of the signal field can be measured, 
from which the quantum state can be obtained in terms 
of phase-space functions [[l0| and the density matrix in 
the Fock basis |nj. The method also applies to matter 
systems and was used to reconstruct the quantum state 
of the (harmonic) center-of mass motion of trapped ions 
]l2t . In this case the coherent displacements are induced 
by rf fields that are resonant to the motional frequency. 

The methods considered so far are restricted to un- 
damped harmonic oscillators. However, in many physi- 
cal systems anharmonic motions (modified by dampings) 
are observed. Recently, interesting nonclassical quantum 
states have been prepared in systems that are typically 
anharmonic. Examples are the generation of amplitude- 
squeezed states in molecules ||[l3| and of Schrodinger 
cat-like states in atomic Rydberg wave packets |l4j. A 
first attempt has been made to reconstruct the quan- 
tum state of anharmonic molecular vibrations using time- 
resolved fluorescence spectroscopy [[[si . It has been 
shown that the density matrix can be obtained by inver- 
sion of high-dimensional systems of linear equations. An 
approach has been given in [ jl6| , extending the pattern- 
function formalism to more general than harmonic poten- 
tials and reconstructing the density matrix from the time- 
dependent position distribution. However, the quantities 
that are typically measured in molecular spectroscopy 
arc different from position distributions in general. For 
example, the time-gated fluorescence spectrum measured 
in experiments of the type described in Q is determined 
by the Franck-Condon overlap factors of the vibrational 
wave functions in two electronic states. For anharmonic 
vibrations this spectrum cannot be identified with distri- 
butions of the type of position distributions and therefore 
it should be used directly in order to reconstruct the den- 
sity matrix p7[. 
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In any case there have been a number of open ques- 
tions, such as those of the determination of suitable sam- 
pling functions mapping the measured data onto the den- 
sity matrix, the choice of optimum observational times, 
and the inclusion into the scheme of damping effects and 
data smearing. The aim of the present paper is to offer 
possible answers to some of them. For all the systems 
mentioned the general problem to be solved is the inver- 
sion of linear equations that relate the measured quanti- 
ties to the density-matrix elements of the system under 
study. This can be done in a very effective way using 
the least-squares method, which dates back to the eras 
of Legendre and Gauss |1S[ ]. Actually, the method has 
been used in the context of quantum-state reconstruc- 
tion problems, such as the reconstruction of the quan- 
tum state of cavity fields by quantum-state endoscopy 
|^9| , vibrations of trapped ion s E^ ], and optical field by 
balanced |2^] and unbalanced Jll[ homodyning. 

For comparison with previous work ||l6| , pif , we will re- 
strict attention to the reconstruction of the density ma- 
trix from the time-dependent position distribution of a 
particle moving in an anharmonic potential. We show 
that the least-squares method can advantageously be 
used in order to reduce the statistical error. We fur- 
ther show that the flexibility of the method enables us 
to perform the reconstruction from data recorded during 
shorter time intervals and to take into account typical ex- 
perimental features, such as smeared and imperfect data. 
Finally, the method can also be used to reconstruct quan- 
tum states of damped systems. It is worth noting that 
the applicability of the method only requires that there 
are measurable quantities that are linear combinations of 
all the relevant density-matrix elements of the quantum 
state of the system under study. 

The paper is organized as follows. In Sec. [n] the prob- 
lem of construction of sampling functions is considered. 
Sections III and [V , respectively, are devoted to the ques- 
tions of observational time and data smearing. Conclud- 
ing remarks are given in Sec. [v|. Elements of the least- 
squares method and the statistical-error analysis are out- 
lined in Apps. |a| and |b[ 



II. SAMPLING FUNCTIONS 

Let us consider a quantum-mechanical system and as- 
sume that at some initial time it is prepared in a state 
with density matrix g n , n ' = {n\g\n'), where \n) are the 
energy eigenstates of the system Hamiltonian. Further, 
let us assume that there is a measurable time-dependent 
(probability) distribution p(x, t) of a quantity x that can 
be given by a linear combination of all density-matrix 
elements g n , n ' that are initially excited, with linearly in- 
dependent coefficient functions S n ^ n i(x,t), 



p(x,t) 



(i) 



When we consider, e.g., a particle that moves in a po- 
tential well and is initially prepared in a bound state 
(e.g., a molecular vibration in a vibronic system below 
the dissociation level), only the discrete part of the en- 
ergy spectrum is excited (n = 1,2,3,.. .). For the sake 
of transparency, in what follows we restrict attention to 
discrete spectra. However, replacing in Eq. (]l|) the sums 
with integrals (or combinations of sums and integrals), 
excitations of continuous parts of the spectrum can be 
treated accordingly. For any physical state the density- 
matrix elements Q n ,n' must eventually decrease indefi- 
nitely with increasing n(n'). Therefore it follows that 
the expression on the right-hand side of Eq. (|IJ) can al- 
ways be approximated to any desired degree of accuracy 
by setting 



Q n ,n> ~ for ra(n') > n E 



(2) 



if "max is suitably large. From the assumptions made 
it is is clear, that Eq. (|l|) can, in principle, be inverted 
in order to obtain the quantum state from the measured 
function p(x, i). 

A powerful method for solving the problem is least- 
squares inversion. To illustrate the method, let us us sup- 
pose that p(x, t) is a (one-dimensional) position probabil- 
ity of a particle moving in a potential well. When during 
the time interval of observation the temporal evolution of 
the quantum state of the particle is only governed by the 
system Hamiltonian, then Eq. ([!]) together with Eq. (|^) 
can be written as 



p(x,t) 



n.n' <n n 



Sn,7i' ^) Qn,n' j 



(3) 



where S n , n '( x > t) is given by 

S n , n >(x,t) = i> n {x)i> n ,{x)e- i ^-^ t . (4) 

Here u> n —cj n i are the transition frequencies of the system, 
the energy eigenfunctions in the position representation 
being given by ip n (x). 

A. Harmonic oscillators 

The simplest system is an undamped harmonic oscil- 
lator of frequency lu^. The equidistant energy spectrum, 
w n — u) n > = (n — n')wh, enables us to separate single di- 
agonals of the density matrix. Introducing the Fourier 
transform 



(5) 



where ojk = kuj\ 1 , k=0, 1, 2, . . ., and T=2tt, we can express 
p( k '(x) in terms of Q n ,n+k as 



p {k \x)= ^ 4'n+k(x)^n(x) Q n +k,n 
n=0 



(0) 
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[^ n (a;) = (v / 7i : 2"n!) _1 / 2 exp(— i 2 /2)H n (i), x being dimen- 
sionless]. Using the least-squares method, Eq. (^|) can 
easily be inverted to obtain the density-matrix elements 
in terms o f th e measured quantit ies. Comparing Eq. (|^) 
with Eq. (Al) and applying Eq. (A7), the reconstructed 
density matrix elements g n +k,n are given by 



Qn+k,n= / dxK^(x)p {k) {x), 



where 



(7) 



(8) 



Here, the matrix F^ fc ^ is the inverse of the matrix G'% 
F( fe ) = (G^) -1 , with the matrix G^ fc ^ being defined by 



G« = / dxBg>{x)Bto(x) 



(m,n = 0, 1, . . . n max -fc), and 

B[ k \x) = ip l+k (x)ipi(x). 



From Eqs. flgj) - (|10|) it follows that 

dxK^{x)B^){x)=5 n>n ,. 



(9) 



(10) 



(11) 



It was found p0| that when n <C fi max then the inte- 
gral kernels (|J) are essentially identical to the sampling 
(pattern) functions derived in B. Moreover, it can be 
shown that only for values of n(rr) that are close to n, max 
the two methods yield substantially different sampling 
functions, the oscillating behavior of the least-squares 
sampling functions being less pronounced than those in 
H . This suggests that the statistical error of the recon- 
structed density-matrix elements Q n , n i with n(n') close to 
"max is smaller for the least-squares method than for the 
method in |J] , because the statistical fluctuation depends 
on the squares of the sampling functions [see Eq. (B4)]. 
The decrease of the statistical error can be understood as 
a consequence of the a priori information on the state to 
be reconstructed: the reconstructed elements Q n ,n' with 
n(n') close to n max cannot be "contaminated" by neigh- 
boring elements with indices above n max . Clearly, when 
the probability of finding excited levels above n max is 
not negligibly small, then the least-squares method can 
cause a systematical error. Taking into account the com- 
plete set of density-matrix elements, we have, in place of 
Eq. (0), 



p {k) (x) = ip n+k (x)ip n (x) g n+k , n 

n=0 

= E + E \ i>n+k(x)ll) n {x) Q n +k,n- (12) 



Using Eqs. (0) and @, we derive 

Qn,n-\-k — Qn,n-\-k 

oo „ 

n' — n max — 

(13) 

where the second term represents the systematical error. 

B. Anharmonic systems 

In order to treat the more general case of non- 
equidistant energy levels, Eq. (^J) can be generalized to 



p( k) (x) = 



n' Qn,n' ) 



(14) 



where p^ (x) is defined by 



p {k) {x) = lim - / dt 



p(x,t). 



(15) 



We see that it is impossible in general to separate sin- 
gle diagonals of the density matrix by integrating the 
position probability over some finite time interval. Actu- 
ally, integration over an infinite time interval needs doing 
in order to exactly single out terms oscillating at chosen 
transition- frequencies of the system, whic h is of course il- 
lusory. We will study this problem in Sec. Ill in more de- 
tail. Here we assume that the terms are (approximately) 
singled out and the density-matrix reconstruction can be 
performed inverting Eq. ([14]) in the same way as for the 
harmonic oscillator, i.e., 



Qn,n> = I dxK n ^(x)p ( - k \x), 



where 



E 



F 1 

± n.n' 



' ^rn , m ' (^) ■ 



(16) 



(17) 



The double-indices n,n' (to, to') are to be chosen such 
that w„( m) - w„'( m ') = wj„, and n, n' (to, to') < n max . For 
chosen k the matrix F is given by F = G _1 , where 



G 



n.n' :m.m' 



B n ,n'{x) = 1p n (x)lp m (x) 



(18) 



(19) 



To give an example, let us consider the bound motion 
of a particle in an anharmonic potential of the Morse 
type, 



n— n—n n 



= ± (e- 



1 



(20) 
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(a > 0). There are Um bound states, where nu = 
[[a~ 2 — 1/2]], with [[y]] being the integral part of y. Their 
wave functions are given by i() n (x)=N n e~ z / 2 z b ' 2 L^ l (z) (n 
= 0, l,...n M ), where z = 2a~ 2 e^ ax , 6 = 2a~ 2 -2n-l, N 2 
= abn\/T(n + 6+1) p2| . Restricting attention to bound 
states we can work with n max <tim- In Fig. [I] we have 
plotted examples of sampling functions K n (x) = K n ^ n (x) 
that are required for the determination of the diagonal 
density-matrix elements g nn [i.e., ui^ = in Eqs. (Ill) 



and (|15|)1, including into the calculation all bound states 
(i.e., rt max = um), which exactly corresponds to the case 
considered in pi[ applying the irregular wave-function 
method (IWM) given in p6| . Comparing the sampling 
functions with those of IWM, from Fig. [jjwe see that the 
latter yields sampling functions that show substantially 
larger oscillations than those obtained by means of least- 
squares inversion. In particular, we see that the effect 
is much stronger than for harmonic oscillators. A re- 
construction based on least-squares inversion is therefore 
expected to give rise to substantially smaller statistical 
fluctuations than IWM. 

To demonstrate this, we have performed computer sim- 
ulations of measurements and reconstructed the diagonal 
density-matrix elements from a set of 5 x 10 3 measured 
data, on assuming the system is initially prepared in a 
state (n\ip) oc a™(n!) _1 / 2 . The exact density-matrix ele- 
ments g nn and the exact time-averaged position distribu- 
tion are shown in Figs. ||(a) and §(b), respectively. The 
position distribution reveals a structure of two peaks lo- 
cated at the turning points. The peak on the side of 
weaker potential is broader than the peak in the region 
of strong repulsive force. Comparing Figs. ||(c) and g(d), 
we see that the reconstruction based on least-squares in- 
version [Fig. ||(c)] yields much less fluctuating results (es- 
pecially for larger values of the quantum number n) than 
that using IWM [Fig. |(d)] . The reason is that in the 
relevant x interval, in which the position distribution is 
essentially nonzero, the least-squares sampling functions 
become weakly oscillating around zero when the quantum 
number n is increased, whereas larger values are attained 
at such x values for which the position probability is small 
[cf. Figs. |l| and ||(b)]. Hence, the least-squares method is 
suited for extracting the information about the density- 
matrix elements from the position distribution even when 
the position distribution is not measured exactly. This is 
not the case when the IWM is used, since the associated 
pattern functions rapidly oscillate with large amplitudes 
over the whole region of probable positions, so that the 
position distribution must be measured with high preci- 
sion in order to extract from it the relevant information 
on the density-matrix elements. 



III. OBSERVATIONAL TIME 

Let us now turn to the problem of measurement time. 
In the case of an undamped harmonic oscillator the sit- 



uation is simple. Since the motion is periodic, observa- 
tion of the position distribution over one period T must 
yield all information about the quantum state. More- 
over, taking into account the symmetry p(x,t + T/2) = 
p(—x, t), we see that observation over one half period is 
already sufficient for the quantum-state reconstruction, 
so that T/2 can be chosen as observational time. In gen- 
eral the excitation energies are not equidistant, and they 
are not discrete even for a bound system. Anharmonici- 
ties prevent the energy levels from being equidistant and 
dampings that accompany any motion give rise to line 
broadenings. 

In order to approximately apply IWM to the density- 
matrix reconstruction of an anharmonic bound system, 
it has been suggested to choose the observational time 
T in Eq. (JsJ) such that it is long compared with all in- 
verse transition frequencies fl6|] . For a Morse oscillator 
observational times that correspond to fractional revivals 
have been proposed To be more specific, the first 

fractional revival for a Morse oscillator appears for T\ 
= 2ir(n M + l/2)/0, where fi = (1 - a 2 /2). Since the 
proposed observational times can be comparable or even 
longer than the characteristic damping times, terms os- 
cillating at discrete frequencies could not be singled out 
in this way. The question raises as to whether or not it 
is possible to reconstruct the density matrix from data 
measured during a substantially shorter time interval. 



A. Factorable sampling functions 

To answer the question, we recall that owing to the 
finite value of n max there is only a finite number of ex- 
ponential functions 



9k(t) 



■ — ^cUfc £ 



-<*v) 



(21) 



that - as long as dampings can be neglected - govern the 
time evolution of p(x,t) [cf. Eqs. @ and (Q)]. We can 
then construct sets of functions that are biorthonormal 
to these exponentials in a chosen time interval in order 
to appropriately decompose p(x,t) and extract from the 
decomposition the density-matrix elements. In principle, 
any interval (small compared with the damping times) 
can be used. Moreover, there are various possibilities 
of constructing a biorthonormal system to a finite set 
of linearly independent functions on a given interval. A 
possible way is to construct them as linear combinations 
of the exponentials gk {t) , where the expansion coefficients 
can be calculate d u sing le ast- squares inversion. From 
App. |A| [see Eqs.Q and @] in a very similar way as 
in the preceding section the orthonormal set of functions 
can be given by 



(22) 



where F = G 1 and the matrix G reads as 
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G k ,i = / dtg* k (t) gi (t), 



(23) 



with T being the length of the chosen time interval. From 
the construction of the functions fk (t) it is clear that 



dt fk(t)9k'{t) = h,k'- 



(24) 



(25) 



In place of Eq. (Ua) we have 



pW(a;)= / dt f k (t)p(x,t), 
Jo 

and the reconstructed density matrix is then 



Qn,n' = I dx J dtK n , n '(x)fk(t)p(x,t), (26) 



with K n>n ' (x) from Eq. (nj 



B. Nonfactorable sampling functions 

The matrix G in Eq. ( f23| ) depends on the chosen time 
interval. When the time interval is too short, then the 
functions gk(t) tend to linearly dependent functions. The 
functions fk (t) become strongly oscillating with large am- 
plitudes, so that even small experimental inaccuracies 
can give rise to big errors. To get reasonable values of 
the statistical fluctuations, the required interval of time 
integration may be too large. The observational time 
however can be drastically reduced when the inversion of 
Eq. (^|) is performed in time and space simultaneously. 
Direct application of least-squares inversion to Eq. (||) 
yields 



Q n ,n> = I dx I dtK n>n >(x,t)p(x,t), 



(27) 



where the time- and position-dependent (nonfactorable) 
sampling functions are given by 



Kn,n' (x, i) 



E 

rn,m'<n n 



Fn,n' ;m.m' ^rn m' > ^) (^8) 



and F = G 1 , with the matrix G being now defined by 

G m ,m'-,n,n' = j ^ X J dt S m m / (x , t) S n , n ' (x , t) , (29) 

where the functions S ntn >(x,t) are given by Eq. (|J). 

Results of reconstruction of the density matrix within 
the framework of Eqs. ( p7| ) - (29|) are plotted in Fig. || 
for the same system as in Fig. |2. In our computer exper- 
iments we have assumed that measurements at N t — 120 
times in a (relatively small) time interval T = 6it / (uji—luq) 
are performed and JV c = 5x 10 3 events at each time are 
recorded. The grid of measurement points is chosen such 



that the systematical error owing to discretization is re- 
duced below the statistical one. It should be noted that 
for chosen value of N t both the systematical and the 
statistical errors of the off-diagonal density-matrix ele- 
ments increase with the "distance" from the diagonal 
elements [for the statistical error, see Fig. |(c)]. From 
Figs. ||(a) -|3](d) we see that least-squares reconstruction 
yields a good estimation of the density-matrix elements 
even for measurement times much shorter than the (first) 
fractional-revival time « 247r/(wi — u>o). In particular, 
Fig ||(d) reveals that the accuracy of the reconstruction 
is comparable with the accuracy that can be achieved - 
for a comparable number of total events - in the case 
when the observational time is extended to infinity. 

It should be mentioned that the reconstruction scheme 
may advantageously be applied also to harmonic oscil- 
lators. An example may be balanced homodyne detec- 
tion of radiation-field modes. Here the phase interval in 
which the quadrature-component distributions are mea- 
sured can be reduced below a ir interval. Recently it was 
suggested to transfer the formalism of density-matrix re- 
construction of quantum harmonic oscillators to classi- 
cal tomographical reconstruction of objects with space- 
varying transparencies p3fl . This can be advantageous 
when the probing beam is very weak and a relatively 
small number of data is available. Application of the 
least-squares inversion with shorter intervals of the an- 
gular variable can make the method also suitable for to- 
mographical reconstruction of objects whose projections 
on some directions are not available. 

It is worth noting that the assumption that the mea- 
surements are performed over the whole x axis may also 
be dropped. When the available x interval is limited, 
then - in close analogy to limited time intervals - the x 
integrals in the reconstruction formulas can be reduced to 
this interval. Hence, the least-squares inversion method 
enables one to reconstruct the density matrix also in cases 
when the time and "position" intervals are smaller than 
the theoretically allowed ones. 



C. Inclusion of damping 

So far we have assumed that damping can be disre- 
garded on the chosen time scale. However, the method 
can also be extended to damped systems. In this case the 
dependence on time of the coefficient functions S n-n i(x, t) 
in Eq. ([!]) is given by appropriately decaying functions in 
place of purely oscillating exponentials exp[— i{oj n — w n i)t\. 
Suppose that the density matrix evolves according to 
some master equation 



Q = £q, 

where £ is a linear superoperator, 



in 



Tig, 



(30) 



(31) 
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with H describing the effect of (Markovian) damping. 
Since the solution of a master equation (|3^) can always 
be represented in the form of 



and 



'(<) 



u, 



m,m \n,n 



(t) On 



(32) 



g n , n i = Qn.n'(O), the above given reconstruction formulas, 
such as Eqs. ( |27| ) - ( p9|) remain valid when the functions 
(x, £) given in Eq. (0) are replaced with 



E 



l/j m (x)l/} m r(x)Ui 



m,m ;n.n 



(t). (33) 



It should be pointed out that in some cases it may 
happen that the G matrix is quasi-singular so that it 
cannot be inverted in the usual way. Physically it means 
that the available data do not carry enough information 
about some of the density-matrix elements. In this case 
regularized inversion can be used, as discussed in App. |a] 
[see Eqs. (|A9|) and (|A10|)]. 



IV. DATA SMEARING 

In general, an x measurement can only be performed 
with finite accuracy and the time cannot be controlled 
precisely, so that in practice one always deals with more 
or less smeared data. For example, in optical homo- 
dyne tomography nonperfect detection yields smeared 
quadrature-component distributions. In molecular emis- 
sion tomography time-resolved fluorescence spectra 
are measured, which necessarily implies that a perfect 
resolution of frequency and time is impossible. Whereas 
attempts have been made to compensate for nonper- 
fect detection in optical homodyne tomography [HQ , an 
open problem has been the inclusion of data smearing in 
quantum-state reconstruction for more general systems. 

Mathematically, smearing means that instead of p(x, t) 
a convolution 



is typically measured, where V(t) and W(x) are single- 
peak functions (centered at zero) describingthe time and 
position windows. Recalling Eqs. (||) and (Q), p(x,t) can 
be related to the density-matrix elements as 



p(x,t) = ^ S n<rl '(x,t) 



Qn,n' i 



n,n'<n n 



where 



with 



S n ,n'(x,t) = V n<n '(t)W n!n '(x), 



V n , n >(t)= / dt'e-^-^'vit'-t) 



(35) 



(36) 



(37) 



W n , n .{x) = / dx'Mx')i>n>{x')W{x' - x). (38) 



dx I dtK nin r(x,t)p(x,t), 



The inversion of Eq. (p5|) can be done in the same way 
leading to Eqs. @-(E9^ i.e., 



(39) 



where K ni7l i (x, t) is calculated according to Eq. ( |28| ) [to- 
gether with Eq. (|29l)l, but with S nt n'(x,t) in place of 
S n<n >(x,t). The limits of integration in Eq. ( p"9| ) are 
given by the chosen intervals of measurement. In practice 
p(x, t) is usually measured on a grid of points {xi, tk}, so 
that the integrals over x and t in Eq. (|39| ) (and the inte- 
grals determining the G matrix) are replaced with sums. 

It should be noted that due to data smearing it may 
be necessary to perform a regularized inversion as men- 
tioned above. For example, when the measured data are 
relatively insensitive to temporal changes of the quantum 
state, then it may be better to set some reconstructed 
density-matrix elements close to zero rather than to let 
them become very large, "trying to fit the data" as close 
as possible. Clearly, one must correctly interpret the re- 
sult - instead of claiming that the elements are measured 
to be zero one should say that there is not enough evi- 
dence for nonzero values. 

The application of the method to systems with data 
smearing is illustrated in Fig. ^ for the same system as in 
Fig. |j| In the computer simulations of measurements the 
data are assumed to be smeared over Gaussian windows, 
V(t) = exp[-t 2 /(2a%)] and W(x) = exp[-x 2 / (2al)}, with 
cr f = 0.27T / (u>i — luq) and a x = 0.3. It is further assumed 
that during an observational time T = 611/(101 — u>o) mea- 
surements are performed at N t = 30 (equidistant) times 
and iV x = 15 (equidistant) positions in an interval —2 < 
x < 10, the total number of recorded events being TVtot = 
10 5 . Hence, at a given point (xi,tk) of the chosen grid 
the number of recorded events is approximately given by 



p(x,t)= dx' dt'V(t'-t)W(x'-x)p(x',t') (34) n hk 



T 



dx / dt W(x - xi)V(t - t k )p(x, t). (40) 



Since owing to time smearing fast oscillating terms can 
hardly be resolved, the measured distribution p(x,t), 
Eq. |3^), becomes insensitive to off-diagonal density- 
matrix elements that are far from the diagonal. To ob- 
tain reasonable results, regularized inversion of Eq. (|35| ) 
is necessary. Using the Tikhonov regularization (see 
App. ^), the regularization parameter used in Fig. ||is A 
= 2 x 10~ 3 . It should be noted that with increasing value 
of A the statistical error [Fig. |J(b)] decreases, but the 
systematical error [Fig. f|(c)] increases. Suitable values 
of A can be obtained from the L curve (see App. |a|) of 
the problem which is plotted in Fig. |[ The best values of 
A correspond to points near the corner. Using singular- 
value decomposition in place of Tikhonov regularization, 
similar results can be obtained. 
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A comparison with results based on measurements on 
the whole time and position scales and without data 
smearing is shown in Fig.^(d) for the diagonal density- 
matrix elements. We see that for smaller values of n the 
results are almost equally good, however with increas- 
ing n the statistical error of the density-matrix elements 
reconstructed from the smeared data strongly increases. 
This can be understood from the fact that typical fea- 
tures (such as rapid oscillations and larger values of x) of 
higher-excited state contributions to p(x, t) are less avail- 
able from the smeared data confined to a limited interval. 

In the examples given above we have included in the 
reconstruction of the density matrix all bound states of 
the anharmonic oscillator (n max = um)- If there is some 
a priori knowledge on the quantum state of the system, 
it may be possible to choose n max such that n max <Um- 
In this way the dimension of the matrix that must be 
inverted can be reduced. An advantage is that the sta- 
tistical error can also be reduced. To illustrate such a 
case, in Fig. |^ we have applied the reconstruction scheme 
to a Morse oscillator whose anharmonicity (a = 0.15) is 
smaller than that in Figs. |^-|^, so that n^=43. Assum- 
ing that the system is again prepared in a state of the 
form considered in Figs. ^ - ^, with the same parameter 
a = — 1.5, we have set n max = 9. With increasing n max 
the systematical error (owing to truncation) can be de- 
creased, but the statistical error is increased. Therefore, 
for given number of available data (and given state) one 
expects that there is an optimum value of n max for which 
the systematical error is just below the statistical one. 



V. SUMMARY AND CONCLUSIONS 

We have considered the problem of reconstruction of 
the density matrix g n . m from the data available in realis- 
tic experiments. The studied systems can be linear oscil- 
lators as well as more general systems with non equidis- 
tant energy spectra. To obtain the complete quantum 
state, the measured quantities must be related to all 
(nonvanishing) density-matrix elements that must nec- 
essarily be known for characterizing the state. Then the 
problem to be solved is the inversion of sets of linear equa- 
tions. Having enough experimental data, such a system 
of linear equations can be overdetermined with respect 
to the density-matrix elements sought. This enables one 
to perform the inversion in different ways. Here we have 
applied the least-squares method. The density-matrix 
elements are obtained as linear combinations of the ob- 
served quantities, so that they fit the measured data as 
close as possible. The main features of the method can 
be summarized as follows. 

(i) The application of the method is not restricted to 
the reconstruction of the density matrix from certain 
specific quantities, such as the time-dependent position 
distribution considered in IWM. Least-squares inversion 
can always be applied when there are measurable quan- 



tities that can be related to all density-matrix elements 
that significantly contribute to the quantum state of the 
system. Moreover, it is not necessary that the system 
evolves undamped, as it is required for applying IWM. 

(ii) The flexibility of the method enables one to take 
into account typical experimental effects and necessities, 
such as data smearing and restrictions that limit the size 
of measurement intervals and the amount of available 
data. For example, using least-squares inversion and 
reconstructing the density matrix of a particle moving 
in a Morse potential from the time-dependent position 
distribution, the measurement time can be substantially 
shorter than in IWM method. 

(iii) Both the reconstruction of the density matrix and 
the estimation of the statistical error can be performed in 
real time. It is worth noting that least-squares inversion 
can advantageously be used in order to reduce the statis- 
tical error below the level given by IWM for a comparable 
amount of data. 

(v) If the measured data are insensitive to certain 
density-matrix elements, so-called regularized solutions 
can be constructed. That is to say, the reconstructed 
values of those density-matrix elements are set close to 
zero instead of dealing with strongly fluctuating values. 
Rcgularization decreases the statistical error of the re- 
constructed density-matrix elements, but simultaneously 
causes a systematical error. It is therefore necessary to 
optimize the regularization such that the introduced sys- 
tematical error is below the statistical one. 

(vi) The method is essentially based on the fact that 
for any physical quantum state the density-matrix ele- 
ments Qn^m must eventually decrease indefinitely with 
increasing n(m), so that the density matrix can effec- 
tively be truncated at some large (but finite) value of 
n(m). Clearly, truncation always introduces a systemat- 
ical error. In practice it is sufficient to choose a trunca- 
tion parameter for which - similar to regularization - the 
systematical error is smaller than the statistical one. 

Finally, it should be mentioned that there have been 
other approaches to the problem of reconstruction of the 
quantum state from finite sets of measured quantities and 
recorded data p5|-|27[|. In particular, the entropic recon- 
struction scheme studied in 1 25 2(| may be used to sys- 
tematically reconstruct the density matrix, starting from 
only a few number of quantities and data and extending 
the scheme to larger sets. In the method suggested in 
p7| conditions are imposed on the reconstruction proce- 
dure such that the reconstructed object is really a density 
matrix (i.e., the diagonal elements must not be negative 
and the trace is equal to unity). Note that the least- 
squares method - and other methods that are based on 
linear transforms of measured data - can produce "nega- 
tive probabilities" resulting from experimental inaccura- 
cies. The reconstruction schemes in p5|-p7f are based on 
the solution of sets of nonlinear equations whose solution 
may become extremely difficult when large sets of data 
must be processed. In contrast, the least-squares method 
enables one to reconstruct the density-matrix elements in 
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real time in a very straightforward way. Such excesses as 
"negative probabilities" can be kept within the statistical 
error bars whose magnitude can easily be estimated and 
eventually decreased by increasing the number of mea- 
surements. 
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APPENDIX A: ELEMENTS OF LEAST-SQUARES 
INVERSION 

To give a summary of least-squares inversion p8| , p0f , 
let us consider a (possibly unknown) no dimensional 
"state" vector f and a stochastic linear transform 



A f 



(Al) 



yielding an mo dimensional (mo > no) "data" vector y 
available from measurements. Here A is a given mo x no 
matrix, and n is an mo dimensional vector whose ran- 
dom elements with zero means and covariance matrix 
W _1 describe the noise associated with realistic mea- 
surements. The task is to infer f from y. 

From the point of view of Bayesian inference, proba- 
bility distributions for y and f can be introduced, and 
the probability for a state vector f under the condition 
that there is a date vector y can be given by 



P(f|y) cx P(y|f)P(f), 



(A2) 



where P(y|f ) is the corresponding conditional probability 
for the date vector y, and P(f ) is an a priori probability 
for the state vector f . We then look for the state vector f 
for which P(f |y) is maximum. Assuming that the noise 
is Gaussian, we have 



P(y|f) cx exp[-i(y - Af)tW(y - Af)] . 



(A3) 



The a priori probability P(f) represents our knowledge 
of the state when no data are available. Hence P(f ) can 
be set constant if no state is to be preferred. Under the 
conditions made maximization of P(f |y) is equivalent to 
minimization of 



C(f) = (y - Af)tw(y - Af), 



from which we see that the minimum is attained at f 
satisfying the equation 



(A4) 
f 



A f WAf = A'Wy. 



(A5) 



If W is diagonal (i.e., the noise is uncorrelated), then 
Eq. (A4) represents a sum of weighted squares of the 
differences between the components of the data vector y 
and the components of the transformed vector Af , each 
term of the sum being multiplied by a weight given by 
the corresponding element of W. Provided that A^WA 
is not singular, from Eq. (|A5|) we obtain 



f = (A t WA)- 1 A t Wy. 



(A6) 



Otherwise, the inversion of Eq. ( A5) is not unique and 
further criteria must be used to select a solution. When 
the matrix W is not known, the n it may be set a multiple 
of a unity matrix, so that Eq. (A6) reduces to 



f = (AtA)- 1 Aty, 



(A7) 



provided that A T A is not singular. Equation (A7) still 



gives the correct averaged inversion, but the statistical 
fluctuation of the result may be (slightly) enhanced. 

If the data are not sensitive enough to some state- 
vector components, then these components can hardly 
be determined with reasonable accuracy. Mathemati- 
cally, A'WA becomes (quasi-)singular and regulariza- 
tions, such as Tikhonov regularization and singular- 
value decomposition, are required to solve approximately 
Eq. JA5| ). For simplicity let us set W = I, where I is the 
unity matrix. 

Using Tikhonov regularization, it is assumed that in 
the absence of significant data some components of the 
state vector can be preferred by a properly chosen a priori 
probability P(f). Assuming a Gaussian distribution and 
preferring the components close to zero (if information 
about them is lacking), we may write 



P(f) cx exp( 



-iA 2 ftf) 



(A8) 



where the (positive) parameter A is a measure of the 
strength of regularizatio n. M aximi zati on of P(f |y) then 
yields, on recalling Eqs. (A2) and (A3), 



f = (A 2 I + A t A) _l A 1 'y- 



(A9) 



Note that (A 2 I + At A) has only positive eigenvalues and 
is thus always invertible. A possible choice of A is based 
on the so-called L curve, which is a log-log plot of ||f|| 
versus ||Ay||, Ay = y — Af, for different values of A. The 
points on the horizontal branch correspond to large noise, 
whereas the points on the vertical branch correspond to 
large data misfit. Optimum choice of A corresponds to 
points near the corner of the L curve. 

Applying singular-value decomposition, the inversion 
of the matrix (A^A) is performed such that their eigen- 
values whose absolute values are smaller than the (posi- 
tive) parameter of regularization <to are treated as zeros, 
but the inversions are set zero (instead to infinity). This 
operation is called "pseudoinverse" of a matrix, 



f = Pseudoinverse(A 1 'A; (7o)A t y- 



(A10) 
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For (To close to zero the result of Eq. (A10) is similar to 
that of Eq. (A7). With increasing <7o smaller absolute 
values of components of f are preferred. 

The effect of the regularization parameters A and ctq is 
similar. The statistical error of the reconstructed state 
vector f is decreased, but simultaneously bias towards 
zero is produced. Hence, optimum parameters are those 
for which the bias is just below the statistical fluctuation. 
The bias can be estimated, e.g., by Monte Carlo generat- 
ing new sets of "synthetic" data from the reconstructed 
state. From these sets one can again reconstruct new sets 
of f . The difference between the mean value of the states 
reconstructed from the synthetic data and the originally 
reconstructed state estimates the bias. 



APPENDIX B: CALCULATION OF THE 
STATISTICAL ERROR 

Let us suppose that the probabilities for observing a 
physical quantity, pi(s), can be related to quantities fk 
as 



Pi(s) 



k 



A,k(s)fk, 



(Bl) 



where I refers to the values of the quantity to be observed 
and s is some parameter that can be changed during 
the observation. Measuring pi(s) for all I and s values, 
we may identify the measured values pi(s) with a data 
vector, from which - according to Eqs. (A7) or (AE) or 
( |A1C| ) - the quantities fk can be reconstructed, 



f k = J2K k ,i(s)pl(s) 



(B2) 



The measured probabilities can be given by pi(s) = 
ni(s)/N(s), where (for chosen s) N(s) and n;(s), respec- 
tively, are the total number of events and the number 
of events yielding the result I. Assuming that ni(s) has 
approximately a Poissonian distribution with mean and 
variance equal to N(s)pi(s), the mean and variance of fk 
can easily be calculated, namely, 



(Re/ fc ) =^Re[X M ( 



l,s 



and 



(Ms)) 

N(s) 

^Re[Kk,i(s)}pi{s) 



Var(Re/ fc )=5:Re[^( S ) l2VarN(s)] 



(B3) 



l.s 



J2MKkAs)} 



N 2 (s) 

2 Pl( s ) 

N(a)' 



and for the imaginary part accordingly. In practice, the 
unknown exact probabilities pi (s) are replaced with the 
estimates pi(s). 
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FIG. 1. Sampling functions for reconstructing the diagonal density-matrix elements Q n<n of a Morse oscillator (a — 0.279) 
from the position distribution for n — 2 (a) and n = 11 (b); full line: least-squares method; dashed line: IWM. 
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FIG. 2. Reconstruction of the diagonal density-matrix elements of a Morse oscillator (a = 0.279), which is assumed to be 
prepared in a state (n|V>) occ^n!" 1 / 2 , a = -l.5, from N c — 5 x 10 3 recorded events in a simulated position measurement; exact 
density-matrix elements g n ,n (a), exact time-averaged position distribution (b), density- matrix elements reconstructed 
using least-squares inversion (c) and IWM (d). The error bars indicate the predicted statistical error (half error bar corresponds 
to a single standard deviation) . 
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FIG. 3. Reconstruction of the density-matrix elements of the same system as in Fig. g from the data of a simulated position 
measurement, in which at each of N t = 120 equidistant times N e = 5 x 10 3 events are recorded, the overall time interval being 
T = 6n/(uji — cjq)', exact density-matrix elements g n ,m (a), reconstructed (real parts of the) density-matrix elements {5 n , m (b), 
predicted statistical error 5g„ trn (c), comparison of the diagonal elements obtained from the data recorded during the actual 
measurement time T (full lines) with the diagonal elements obtained (for the same total number of events) in the case when T 
— > <x (dashed lines) (d). 
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FIG. 4. Reconstruction of the density-matrix elements of the same system as in Fig. g from a simulated measurement of 
the smeared distribution p(x,t) [at = 0.2 7r/(a>i — u>o), &x = 0.3] at Nt = 30 equidistant times [during the observational time 
T = 67r/(aji — loq)] and JV X = 15 equidistant positions [in the interval —2 < x < 10] for a total number of events of iVtot = 10 , 
using Tikhonov regularization with \ — 2 x 10 -3 ; reconstructed (real parts of the) density-matrix elements § n ,m ( a )> predicted 
statistical error &Q n , m (b), systematical error 5 B Qn,m due to regularization (c), comparison of the diagonal elements obtained 
from the smeared data recorded during the actual measurement time T (full lines) with the diagonal elements obtained (for 
the same total number of events) without data smearing, and T^oo (dashed lines) (d). 
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FIG. 5. The L curve for the problem considered in Fig. |; A = 1(T 4 (o), 2 x 10~ 3 (+), 5 x 10~ 3 (*), 5 x 1(T 2 (x). 
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FIG. 6. Reconstruction of the density-matrix elements of a system of the type considered in Fig. g, but with smaller 
anharmonicity, a — 0.15, and a truncation parameter smaller than the number of bound states, n max = 9, from a simulated 
measurement of the smeared distribution p(x,t) [at = 0.2 n/(u)i — uio), Ox = 0.3] at N t = 30 equidistant times [during the 
observational time T = Att / {uj\ —ujq)] and JV X = 15 equidistant positions [in the interval — 2 < x < 10] for a total number of events 
of JVtot = 10 6 , using Tikhonov regularization with A= 10 -3 ; reconstructed (real parts of the) density-matrix elements g n , m (a), 
predicted statistical error Sg n , m (b). 
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